Vortex states in mesoscopic superconducting squares: Formation of vortex shells 
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We analyze theoretically and experimentally vortex configurations in mesoscopic superconducting 
squares. Our theoretical approach is based on the analytical solution of the London equation using 
Green's-function method. The potential-energy landscape found for each vortex configuration is 
then used in Langevin-type molecular-dynamics simulations to obtain stable vortex configurations. 
Met ast able states and transitions between them and the ground state are analyzed. We present our 
results of the first direct visualization of vortex patterns in /xm-sized Nb squares, using the Bitter 
decoration technique. We show that the filling rules for vortices in squares with increasing applied 
magnetic field can be formulated, although in a different manner than in disks, in terms of formation 
of vortex "shells" . 

PACS numbers: 74.25.Qt,74.25.Ha,74.78.Na 



I. INTRODUCTION 



The growing interest in studying vortex matter 
in mesoscopic and nano-patterned superconductors is 
closely related to recent progress in nano-fabrication and 
perspectives of their use in nano-devices manipulating 
single flux quanta. As distinct from bulk supercon- 
ductors, vortex states in nano- and mesoscopic samples 
are determined by the interplay between the intervor- 
tex interaction (which is modified due to the presence 
of boundaries) and the confinement. In general, the 
shape of a mesoscopic sample is incommensurate with 
the triangular Abrikosov lattice, and as a consequence, 
the resulting vortex patterns display strong features of 
the sample shape and may differ strongly from a tri- 
angular lattice. Strong finite size effects in conjunction 
with strong shape effects determine the vortex configu- 
rations. For example, in mesoscopic disks vortices, as 
shown theoreticall}*iii^i^i^i^i^i^ and experimentally^, form 
circular symmetric shells (similar to two-dimensional 
(2D) system of charged classical particles^). Moreover, 
due to strong confinement effects in small disks vor- 
tices can even merge into a giant vortex (GV), i.e., a 
single vortex containing more than one flux quantump^, 
as was recently confirmed experimentally^^. Further- 
more, it was recently demonstrated^^ that vortices can 
merge into a cluster or a GV in /im-sized mesoscopic 
niobium disks which is induced by strong disorder in 
combination with rather weak confinement, while nei- 
ther of these effects alone would lead to a GV/cluster 
formation. Similarly, shape- and symmetry-induced vor- 
tex patterns can be formed in mesoscopic superconduct- 
ing trianglesi^iii^ii^, squaresi^ii^iiiii^, or, in general, in 
symmetric polygons^-'^. However, unlike disks where 



the vortex patterns result from the interplay between 
the discrete symmetry of the (triangular) vortex lattice 
and the cylindrical (Coo) symmetry of the disk, meso- 
scopic polygons have discrete symmetry that can coincide 
(triangles, Cs symmetry) or include as a subgroup (e.g., 
hexagons with Cq symmetry) the symmetry of the vortex 
lattice. In such cases highly stable vortex configurations 
are possible for some values of magnetic field (providing 
commensurate numbers of vortices) because the vortex- 
vortex interaction is enhanced by the effect of bound- 
aries. Strikingly, strong boundary effects can even lead to 
symmetry-induced vortex states with antivorticesi^Ai^ii^ 
(i.e., the symmetry of the vortex configuration with an- 
tivortices can be restored by the generation of a vortex- 
antivortex pair). 

In contrast to Csn-symmetric (where n is an integer) 
polygons, squares are incommensurate with triangular 
vortex lattice for any applied magnetic field. The vortex- 
vortex interaction and the effect of boundaries are al- 
ways competing in mesoscopic squares. Resulting from 
this interplay: i) the ground state of the vortex system 
always involves nonzero elastic energy and, as a conse- 
quence, ii) there are metastable states with energies close 
to the ground state (or, in principle, the ground state 
even could be degenerate). Early studies on vortices in 
mesoscopic squares were either limited to very small sam- 
ples with characteristic sizes of the order of ^ (where ^ is 
the coherence length) which were able to accommodate 
only few vortices^^, or they focused on the possibility of 
generation and stability of vortex-antivortex patterns in 
squares^^aili^. Here we present a systematic theoretical 
analysis of vortex configurations in mesoscopic squares 
and their first direct observation in /im-sized niobium 
squares using the Bitter decoration technique. To study 
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the formation of vortex patterns and transitions between 
the ground and met ast able states, we analyticahy solve 
the London equation using the Green's function method, 
and perform molecular-dynamics simulations. To obtain 
the stable vortex configurations, we analyze the filling of 
squares by vortices with increasing applied magnetic field 
and the formation of vortex "shells", similarly to those 
observed in disks. 

The paper is organized as follows. The theoretical for- 
malism and the solution of the London equation using the 
Green's function method, for a system of L vortices in a 
rectangle sample, are described in Sec. IL In Sec. Ill, we 
discuss the evolution of vortex configurations with mag- 
netic field calculated using the solution of the London 
equation found in Sec. II and the molecular-dynamics 
simulations (Sec. III. A). We formulate the filling rules 
and discuss the formation of vortex shells in mesoscopic 
superconducting squares in Sec. III.B. Metastable states 
and the transitions between them and the ground state 
are analyzed in Sec. III.C. In Sec. IV, we present the 
results of our direct experimental observations of vortex 
patterns in niobium squares using the Bitter decoration 
technique, and compare the calculated patterns with the 
experimentally measured vortex configurations. The con- 
clusions are given in Sec. V. 

II. THEORY: THE LONDON APPROACH 

We consider a strong type-II superconductor (i.e., 
characterized by the Ginzburg-Landau parameter n = 
A/^ ^ 1, where A is the London penetration depth and ^ 
is the coherence length) with rectangular cross section in 
the x-y plane and thickness d in the 2:-direction. Note 
that the London approach is applicable also for weak 
type II superconductors in case of thin-film samples with 
thickness d <C A where the penetration depth is modi- 
fied: A ^ A = A/d^, or in case of low vortex densities 
in rather large mesoscopic samples (i.e., with the lateral 
dimensions a, a > A) where vortices are well separated 
and the order parameter is = 1 everywhere except 
at the vortex cores. The latter case corresponds to our 
experiments with /im-sized niobium squares as described 
below. In our model the external magnetic field H is 
applied normal to the x-y plane, i.e., along the z-axis: 
h = /iz. We also assume that the vortex cores are straight 
lines along the z-direction. Then the local magnetic field 
can be found by solving the London equation: 

L 

- X^V^h ^h = ^ohY, ^(^ - ^i)^ (1) 

where is the fiux quantum and {r^ = (xi^yi)^ i = 
1, . . . , L} are the positions of L vortices. If we also ne- 
glect the distortion of the external magnetic field due to 
the sample, i.e., assume that the value of the magnetic 
field outside the sample near its boundary is equal to 
the applied field, then the boundary conditions for the 







b 















-a/2 





a/2 



FIG. 1: The cross-section of a rectangular superconductor 
with sides a and b. The external magnetic field H is applied 
along the z-axis, and its value is assumed to be constant out- 
side the sample. 

magnetic field are: 

h{±a/2, y) = h{x, 0) = h{x, b) = H. (2) 

The geometry of the problem is shown in Fig. [H The 
Green's function method for solving the London equation 
Eq. ([1]) with the boundary conditions Eq. ([2|) was previ- 
ously used by Sardella et al}^ . However, they limited 
themselves to the special case where one of the sides of 
the rectangle is much larger than the other, i.e., a stripe. 
Such an approximation considerably simplifies the prob- 
lem but the resulting solution missed the generality (the 
symmetry with respect to the permutation x ^ y) and 
thus could not be used in our case of a square: a = b. We 
seek for a solution of Eq. ^ with the boundary condi- 
tions Eq. ([2]) which is valid for a rectangle with arbitrary 
aspect ratio a/b. The Green's function associating with 
the boundary problem defined by Eqs. ([T]) and ([2]) must 
satisfy the following equation: 

- A^V^G + G = (5(x - x')S{y - y'), (3) 

and the boundary conditions: 

G(±a/2, y) = G(x, 0) = G(x, b) = 0. (4) 

Multiplying Eq. ([I]) by G and Eq. by h and subtract 
one from another, we obtain 

-A2(GV'/i- W'G) 

L 

= G^o ^ (5(r - - hS{x - x')5{y - y'). (5) 

i=l 

Integrating Eq. ([5]) over the sample area, we arrive at 

/a/2 f>h 
dx / dy{G\/^h-h\/^G) 
-a/2 Jo 

/a/2 pb L 

dx / d^[G<I>o V(^(r-ri)- 
-a/2 Jo 

hS{x-x')5{y-y')]. (6) 
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Further we use Gauss theorem, 
ro/2 rb 



/a/ z po 
dx / dy{GV^h-hV^G) 
-a/2 Jo 



boundary 



dl{G—- - h—-), 

^ on on^ 



where d/dn is the derivative in the normal direction to 
the boundary, and the boundary conditions Eqs. and 
([2|), and we find the expression for the magnetic field: 



h{x',y') =H 



/a/2 nb 
dx / dyG{x,y,x\y' 
-a/2 Jo 



a/2 

/2 
L 

^o^G{xi,yi,x\y'). 

i=l 



(7) 



Therefore, the problem of finding the solution for the 
local magnetic field is reduced to the determination of 
the Green's function G{x^y^x' ^y'). In order to find a 
solution to Eq. (j3|) with the boundary condition Eq. (|3|), 
we expand the Green's function in a Fourier series, 

, 2 . ,m7Ty . . .miry. , 

G{x, y^x ,y) = - sm(— ^) sm(— x ). 

m=l 

(8) 

Note that the boundary conditions Eq. (jH) are satisfied 
at = 0, 6. Further we substitute this expansion into Eq. 
(|3]) and obtain 



m=l 



d ym{x^x') . mny' . mny 



dx'^ 



■ sin( 



sin(- 



— ) ^m(^,^ )sm(— ^)sm(— ^) 



. .miry . miry 
+ sm(— ^) sm(— X ) 



5{x -X)- sm(— ^) sm(— ^), 



(9) 



where we used the following J-function representation 
^iy-y) = lZ^ sm(— ^) sm(— ^) 

171=1 

since | Y^sin(^^^^), m = 1,2,3...| forms a complete 

set of orthonormal functions. As a result, we obtain 
the following equation for the Fourier-transform of the 
Green's function 



d'^gmjx.x') 
dx'^ 



+ a^^m(^, x') = S{x - x'), (10) 



where 
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(11) 



The functions gm{x^x') must satisfy the boundary con- 
ditions ^^(±a/2,x') = 0. In order to solve Eq. ([TO]), we 
first take its Fourier transform. 



where 



from which we obtain a particular solution to Eq. ([TO]) 



Qra I a- 



-Oim\x — x'\/X 



[cosh(a^(x — x')/\) — sinh(a^|x — x'\\)\ . 



2a^A 

The general solution of Eq. (p!0|) reads as: 
1 



■ [cosh(am(^ — x')/\) — sinh(Qf^|x — x'\\)\ 



2arnA 

+ A{x') sinh(am^/A) + B{x') cosh(am^/A) 

= — ^— - [ — sinh(am|2; — x'\\) + C{x') sinh(am^/A) 
2amA 

+ D{x') cosh(am^/A)] . 

Using the boundary conditions Eq. (|4]) we find the coef- 
ficients C{x') and D{x'): 

C{x') = — coth(a^a/2A) sinh(x'); 
D{x') = tanh(arn<^/2A) cosh(x'). 

Then the solution for gm{x^x') is given by 

1 



gm {x^X ^ 



2\am sinh(amtt/A) 

X I cosh [a^nd^ — x'\ — a) /a] 

- cosh [am{x + x')/\] \. (12) 



Inserting this result into Eq. (|8]), we obtain the expres- 
sion for the Green's function: 



m=l 



1 



2Aam sinh(Qfrn<^/A) 
cosh \am{x^x')/\\ >. 



cosh [a^(|x — x\ — a)/ a] 
(13) 
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From it we obtain an expression for the local magnetic 
field: 



4 v-^ h 
_| \ i 

^ cosli(a2m+i^/A) 
cosli(a2m+i<^/2A) 



sli[(y - h/2)/\] 
cosli(V2A) 

(2m + l)7ry 



(14) 



Note that this solution is valid for a rectangle with ar- 
bitrary aspect ratio a/h and is a generalization of the 
earlier result presented in Ref. [l^. 



Using the obtained solution or the London equation for the local distribution of the magnetic field h{x^ y), we obtain 
the Gibbs free energy per unit length of an arbitrary vortex configuration: 

^ = ^ /"^./^ie^d ^ A ^ ^core ^ ^field 

i=l ^ j=l ^ 



qH f cosh[{yi - b/2)/\] 4 _2 b . r (2m + l)7rt/i -, cosh(Q2m+ia:i/A) | 

TTyl^t cosh(6/2A) 6^p"2m+i(2TO+ l)7r'^''^ L b ^ cosh(a2„»+ia/2A) J 

_^^^^( f tanh(6/2A) 8 ^ tanh(Q2m+ia/2A) ] 

SttA ^ """"^ ^TT 1 5/2A 7r2 [(2m+ l)a2„.+i]2(a2™+ia/2A) / 
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(15) 



Here, A = a x b is the area of the rectangle. The last 
two terms are the energies associated with the external 
magnetic field and the vortex cores, respectively. The 
Green's function in the first term describes the interac- 
tion between vortices and also the interaction between 
vortices and their images which are situated outside the 
sample. The second term represents the interaction be- 
tween the ith vortex and the shielding currents. Note 
that in Ref. the authors limited their consideration 
to the case of a thin film such that (ttA/^)^ 1 and the 
term "1" in Eq. (pTj) can be neglected. The London the- 
ory has a singularity for the interaction between a vortex 
and its own image (self-interaction) . We notice that when 
i = j the Green's function does not converge. To avoid 
the divergency, we apply a cutoff procedure (see, e.g., 
Refs 2Q,j2lll22l ) . which means a replacement of |r — Vj \ by 
for i = j. It was shown in Ref. [S^ that the results 
of the London theory agree with those of the Ginzburg- 
Landau theory, the vortex size should be chosen as V^^, 
and therefore we take a = \pl. The confinement energy 
is given by tc = ef^*^^'^ + ta. 

In Figs [2]^ a) and (b), we plot the distribution of the 
confinement energy for mesoscopic squares with a = 3A 
and a — 15 A, correspondingly. In the mesoscopic square 
with a — 3A, Fig. [2](a), the screening current extends in- 
side the square and interacts with all the vortices. But 
in the large mesoscopic square (we call it "macroscopic" ) 
with a = 15A, only the vortices which are close to the 



boundary feel the screening current. In the mesoscopic 
square, vortices strongly overlap with each other (see 
Fig. [2]^c)) while in the macroscopic square, the inter- 
action between vortices is rather weak and only the clos- 
est neighbors are important (see Fig. [2^d)). This differ- 
ence between small (mesoscopic) and large (macroscopic) 
squares leads, in general, to the size-dependence of the 
vortex patterns in mesoscopic samples as it was recently 
demosntrated for disks (see Ref. [ol) . 

III. THE EVOLUTION OF VORTEX PATTERNS 
WITH MAGNETIC FIELD 

A. Molecular-Dynamics simulations of vortex 
patterns 

Within the London approach, vortices can be treated 
as point-like "particles", and it is convenient to employ 
Molecular Dynamics (MD) for studying the vortex mo- 
tion driven by external forces (see, e.g., Refs. 0, [ill, HI, 
[25I), similarly to a system of classical particles^. In the 
previous section we obtained the analytic expression for 
the free energy of a system of L vortices as a function 
of the applied magnetic field, Eq. (p!5|) . The force felt by 
the ith vortex can be obtained by taking the derivative 
of the energy: 

F^ = -V,^, (16) 
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where Vi = ^e^^ + ^^2/ dimensional deriva- 

tive operator. 

The overdamped equation of vortex motion can be pre- 
sented in the form: 



■ self 



M ■ 



(17) 



where the first three terms are as follows: Fij is the force 
due to the repulsive vortex- vortex interaction of the ith 
vortex with all other vortices, F^^^j is the interaction 
force with the image, and is the force of interaction 
with the external magnetic field which enters the sample 
through the boundaries; r] is the viscosity, which is set 
here to unity. Note that Eq. ([16]) contains these three 
terms (with the free energy defined by Eq. (p!5|)). and in 
Eq. ([T7j) we added a thermal stochastic term to sim- 
ulate the process of annealing in the experiment. The 
thermal stochastic term should obey the following condi- 
tions: 



and 



{Ff{t))=0 (18) 



It is convenient to express the lengths in units of A, the 
fields in units of i^c2, the energies per unit length in units 
of = ^q/SttA • and the force per unit length in 

units of /o = ^q/SttA-I/ ^ where A is the sample's area. 
In our calculations we use the value of the Ginzburg- 
Landau parameter Hi = 6 taken from the experiment with 
Nb (see below). 

In order to find the ground state vortex configurations 
in squares, we perform stimulated annealing simulations 
by numerically integrating the overdamped equations of 
motion Eq. (pT|) . The procedure is as follows. First we 
generate a random vortex distribution and set a high 
value of temperature. Then we gradually decrease the 
temperature to zero, i.e., simulating the annealing pro- 
cess in real experiments (see, e.g., Ref. '26). To find the 
minimum energy configuration, we perform many simu- 
lation runs with random initial distributions and count 
the statistics of the appearance of different vortex con- 
figurations for each L. This procedure simulates^ the 
statistical analysis of experimental data with simultane- 
ous measurements of vortex configurations in arrays of 
many (up to 300) practically identical samples. It was 
used in experiments with Nb disks in Refs. pTlflll and also 
in experiments with Nb squares presented in this paper. 



FIG. 2: (Color online) The profiles of the confinement energy 
Cc = ef^'^^"^ + en (measured in units of go = ^I/StvA • 1/A^, 
where A is the area of the sample) for mesoscopic supercon- 
ducting squares with size a = 3A (a) and 15A (b). The Gibbs 
free energy distributions for squares with a = 3A (c) and 15 A 
(d) for the vortex state with L = 5. 



B. Filling rules for vortices in squares with 
increasing magnetic field: Formation of vortex shells 

The results for the vortex patterns for different vortici- 
ties L are shown in Figs [3] and [H With increasing applied 
magnetic field, vortex configurations evolve as follows: 
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FIG. 3: (Color online) The evolution of vortex configurations 
for the states with vorticity increasing from L — 1 to 12, 
in a superconducting square with a = 3A (the same results 
found for larger squares, e.g., with a = 15A). The vortices 
in the outer shell are shown by the blue (black) circles while 
the inner-shell vortices are shown by the yellow (grey) circles. 
The formation of the second shell starts when L — 5. 



Starting from a Meissner state with no vortex, the first 
vortex appears in the center - see Fig. [3l^a), for L = 2 
the two are located symmetrically on the diagonal - see 
Fig. [Sjb). Further increase of the magnetic field leads 
to the formation of a triangular vortex pattern having 
a common symmetry axis with the square, which is the 
diagonal - see Fig.[3](c). For L = 4 vortices arrange them- 
selves in a perfect square. Fig. [31(d), whose symmetry is 
commensurate with the sample and therefore it turns out 
that this is a highly stable vortex configuration^^^^. Note 
that even in the bulk the gain in the elastic energy is very 
small during the transition from the triangular vortex lat- 
tice to the square one, and consequently, in the presence 
of a square boundary, it turns out that a square vortex 
lattice can be easily stabilized (for commensurate vortex 
numbers). For vorticity L = 5, vortices tend to form ei- 
ther a pentagon, or a square with one vortex in the center 
(see Fig. [3fe), the transition between this configuration 



and the pentagon-like pattern will be discussed below). 
The additional vortex appears in the center thus form- 
ing a second shell in a similar way as in disks^iii^, but 
in the latter, this occurred for a larger L- value (L = 6). 
To distinguish different shells and indicate the number 
of vortices in each shell, we use the same notations as in 
Refs.'o'.?.'^. For example, the pentagon-like configuration 
and the pattern with four vortices in the outershell and 
one vortex in the center are denoted as (5) and (1,4), 
respectively. (It is clear that vortex shells in squares are 
not as well defined as in disks and sometimes it is a mat- 
ter of choice how to define them.) Compared with disks, 
which have Coo symmetry, the C4 symmetry of squares 
induces a new element of symmetry in the resulting vor- 
tex patterns. In other words, vortex patterns in squares 
(tend to) acquire elements of the C4 symmetry even if 
they are not arranged in a perfect square lattice. For ex- 
ample, the calculated vortex patterns share one (L = 6, 
Fig.[3jf)) or two {L = 7 and 8, Figs.[3jg) and (h), corre- 
spondingly) symmetry axes of the square parallel to its 
side. This tendency to share symmetry elements with 
the square boundary remains also for larger vorticities as 
can be seen, e.g., in Figs.[3Kj), (k)? 0-) vorticities 
L = 10, 11, and 12, respectively. For the commensurate 
number of vortices L = 9, a perfect symmetric square- 
lattice pattern is formed. 

Using the concept of vortex shells, we analyzed the 
filling rules for mesoscopic superconducting squares with 
increasing magnetic field. To summarize these rules, for 
L = 1 to 4, vortices are arranged in a single "shell" ; the 
second shell appears when L = 5, and then vortices fill 
the shells as follows: As the vorticity L increases from 
L = 5 to 9, the new vortices fill the outer shell. Then the 
number of vortices in the inner shell starts to increase for 
L > 9 (see Figs. [3fj), (k), and (1)). This occurs because 
the outer shell is formed by 8 vortices (i.e., three per each 
side) which turns out to be stable. Thus, the new vortices 
fill the inner shell until L = 12. Then, again, the newly 
generated vortices start to fill the outermost shell until 
L = 16, when the number of vortices in the outermost 
shell becomes 12, which is also stable (i.e., commensurate 
with the square boundary). The formation of the third 
shell starts when the vorticity becomes L = 17 (note that 
for L = 17 the vortices can arrange themselves either in 
a two-shell configuration (5,12), or in a three-shell con- 
figuration (1,4,12) which occurs to have a slightly lower 
energy, see analysis below). In a similar way, the filling 
of shells occurs for larger values of L (e.g., for 3-, 4-shell 
patterns, etc.). As a general rule, the outermost shells 
containing AN vortices, where N is an integer, are very 
stable. With increasing the density of vortices, the av- 
erage distance between them decreases. As a result, the 
interaction between vortices becomes more and more im- 
portant leading to the formation of the triangular-lattice 
phase away from the boundary. Therefore, the triangu- 
lar lattice is recovered for large vorticities being distorted 
near the square boundaries. Note that for large enough 
L vortices do not form a square lattice even for commen- 
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FIG. 4: (Color online) The evolution of vortex configurations 
for L = 15 to 18 (a)-(d), and for L = 25 (e) to 29 (f), in a 
superconducting square with a = 3A. For vorticities L = 15 
to 18 ((a)-(d)), the outermost shell formed by 12 vortices 
is complete (commensurate with the square boundary), and 
with increasing magnetic field vortices fill inner shells. Note 
that when the inner shell also becomes complete (L = 16, 
state (4,12) (b)), the third shell starts to form for L = 17 
(c). For states with larger vorticities, e.g., L = 25 (e) L = 29 
(f), the vortex patterns are very close to a triangular lattice 
which is distorted near the boundary. 



surate vortex numbers (e.g., for L = 25, 36, etc.) as it 
does for L = 4, 9, and 16. Some expamples of two- and 
three-shell vortex patterns are shown in Fig. [H 



C. The ground state and met ast able states 

The incommensurability of the square boundary with 
the triangular vortex lattice creates metastable vortex 
configurations. While in many cases metastable states 
are well separated in energy from the ground state, in 
some cases, namely, for borderline configurations having 
n and n+1 shells, the lowest-energy metastable state can 
become almost indistinguishable from the ground state. 
In such cases, vortex states with very close energies can 
have comparable probability to be realized experimen- 
tally. An example of a such state is the case L = 5. The 
stable states for L = 5 are shown in the insets of Fig. [5l 
In order to examine which one is more stable, we inves- 
tigate the free energy as a function of the displacement 
of one of the vortices while we allow the other vortices 
to relax to their lowest-energy positions. We start with 
the pentagon-like configuration (5) (the left inset) and 
we change the position of this vortex moving it towards 
the center of the square and let the other vortices ad- 
just their positions accordingly. At the end, we arrive 
at the square- symmetric state (1,4). We plot the free 
energy of the system as a function of the displacement 
of this vortex from its equilibrium position, and we re- 
peat this procedure for all the vortices A, B, C, D, and 




FIG. 5: (Color online) The change of the free energy (G — 
Go) versus the displacement R of one of the vortices in the 
initial pent agon- shaped configuration from its initial position 
towards the center (two different lines for each configuration 
correspond to increasing and decreasing AR as shown by the 
arrows in (a)). Go is the free energy associated with external 
magnetic field and the vortex cores (term "4" in Eq. ([IS])), 
which is independent of the positions of the vortices. The 
two stable states, the pentagon-like state (5) and the square 
symmetric state (1,4), are shown in the insets. The vortices 
are labeled by A, B, C, D and E. Two different symmetry axes 
of the configuration (5) are shown by the dash-dotted line in 
the insets of (a) and (b), respectively. The side of the square 
is a = 3A. In both cases, the configuration with one vortex 
in the center (1,4) has a lower energy than the pentagon- like 
pattern (5). Note that the curves for B and D (and for A and 
E) are slightly different due to the fact that the configuration 
(5) is not perfectly aligned with respect to the symmetry axes. 



E (we always move only one vortex while all others relax 
to minimize the free energy). For any of the five vortices, 
this procedure leads to a barrier between the two states. 
We notice that there are two possible pentagon-like con- 
figurations (5) which share different symmetry axes with 
the square, see Figs.[5fa) and (b). The difference of their 
free energy is less than 10~^. In Fig. [5]^a) we see that 
the motion of vortex C is accompanied with the lowest 
energy barrier. This is because vortices A, B, D and E 
are already close to their final positions in state (1,4). 
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vortices form the outermost shell and the other five can 
form either a one-shell or two-shell configurations sim- 
ilarly as state L = 5. Again, we move one of the five 
vortices in the inner shell of the state (5,12) to the center 
of the square. The analysis of the free energy shows that 
the difference of the free energy between the two states 
(|AG| ~ 10~^) is much smaller compared to the states 
for L = 5 (|AG| - IQ-^). The reason for this is that for 
1/ = 17, the twelve vortices in the outermost shell can 
adjust themselves to lower the free energy, which create 
much "softer" walls for the five vortices in the inner shell 
than the sample boundary. Thus, the change of the free 
energy due to the movement of the vortices in the inner 
shells can be more or less compensated by the movement 
of the vortices in the outermost shell. 



FIG. 6: (Color online) The change of the free energy (G — Go) 
versus the displacement R of one of the vortex in the inner 
shell of the state (5,12) from its initial position towards the 
center; Go is defined in the caption of Fig. [5] The change 
in the free energy due to the movement of the vortices in 
the inner shells (i.e., (5,12) (1,4,12)) is damped by the 
movement of the vortices in the outermost shell which act as a 
"softer" wall than the boundary (in the case of the transition 
(5) (1,4), see Fig. (5]). The movement of the vortices in 
the outmost shell causes more saddle points. The two states, 
(5,12) and (1,4,12), have very close free energies. 



Moving vortices B or D lead to a higher energy barrier. 
Finally, moving vortex A or E to the center is associated 
with the highest barrier and passing over a saddle point 
(jump in G — Go)- Then we move the central vortex of 
state (1,4) back to its initial positions in state (5). The 
highest-barrier transitions (i.e., curves A and E) show a 
hysteretic behaviour which is an indication of metastable 
states. 

In Fig. [51(b), we show the results of the calculation of 
the free energy as a function of displacement of a vortex, 
for a different modification of the state (5), i.e., when 
the vortex configuration has the symmetry axis coincid- 
ing with the diagonal of the square (cp. Fig.[5fa)). Note 
that these two configurations of state (5) have practi- 
cally the same free energy and thus equal probability to 
appear in experiment. Moving vortex E, which is situ- 
ated on the diagonal of the square (see the left inset in 
Fig. [SJb)), is accompanied by the highest energy barrier 
compared to moving other vortices. The reverse process 
(i.e., moving the central vortex to position E) leads to a 
very high potential barrier, and the pentagon-like state 
cannot be restored unless a random (thermal) force is 
added to break the symmetry. Moving vortex B or C is 
accompanied by the lowest energy barrier. State (1,4) 
has a lower free energy than state (5). According to our 
calculations, it is the ground state for L = 5. 

Similar transitions are found between two- and three- 
shell vortex configuration for L = 17 (see Fig. [6]). Twelve 



IV. EXPERIMENTAL OBSERVATION OF 
VORTEX CONFIGURATIONS IN MESOSCOPIC 
NB SQUARES 

To visualise the corresponding vortex configurations 
experimentally we used the well-known Bitter decora- 
tion technique which is based on in situ evaporation 
of 10 — 20 nm Fe particles that are attracted to re- 
gions of magnetic field created by individual vortices 
and thus allow their visualisation (details of the tech- 
nique are described elsewhere^^). The mesoscopic sam- 
ples for this study were made from a 150 nm thick Nb 
film deposited on a Si substrate using magnetron sputter- 
ing. The film's superconducting parameters were: transi- 
tion temperature Tc = 9.1 K, magnetic field penetration 
depth A(0) ~ 90 nm; coherence length ^(0) ~ 15 nm; 
upper critical field Hc2{0) ^ 1.5 T. Using e-beam lithog- 
raphy and dry etching with an Ar ion beam through a 
250 nm thick Al mask, the films were made into arrays 
of small square "dots" of 4 different sizes, with the side 
of the square, a, varying from 1 to 5 /im. Each array 
typically contained 150 to 200 such dots. A whole array 
was decorated in each experiment, allowing us to obtain 
a snapshot of up to 100 vortex configurations in dots of 
the same shape and size, produced in identical condi- 
tions (same applied magnetic field H and temperature 
T, same decoration conditions). It was therefore possible 
to simultaneously visualise vortex configurations for sev- 
eral different vorticities L (in samples of different sizes) 
and also to gain enough statistics for quantitative analy- 
sis of the observed vortex states in terms of their stabil- 
ity, sensitivity to sample imperfections, and so on. Be- 
low we present the results obtained after field-cooling to 
T 1.8 K in perpendicular external fields ranging from 

= 20 to 60 Oe. We note that the above temperature 
(1.8 K) represents the starting temperature for the ex- 
periments. Thermal evaporation of Fe particles usually 
leads to a temporary increase in temperature of the dec- 
orated samples but the increase never exceeded 2 K in 
the present experiments, leaving the studied Nb dots in 
the low-temperature limit, T < 0.5 Tc. 
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FIG. 7: Scanning electron microscope (SEM) images of 
vortex configurations observed experimentally for vorticities 
L = 2 to 13. Vortex positions are indicated by clusters of 
small white (Fe) particles, (a) L = 2; sample size (side of 
the square) a ^ 2.5 /xm, H = 20 Oe; (h) L = 3; a ^ 2 /im, 
if = 35 Oe; (c) L = 4; a ^ 2.4 /xm, if = 40 Oe; (d) L = 5; 
a ^ 2.4 /im, if = 40 Oe; (e) L = 6; a ^ 2.5 /xm, if = 40 Oe; 
(f) L = 7; a 2 /im, if = 60 Oe; (g) L = 9; a ^ 3.5 /xm, 
fT = 35 Oe; (h) L = 10; a ^ 3.5 /urn, fT = 35 Oe; (i) L = 10; 
a ^ 3.5 /xm, if = 35 Oe; (j) L = 11; a 2.5 /xm, if = 60 Oe; 
(k) L = 12; 2.6 /xm, fT = 60 Oe; (1) L = 13; a 5 /xm, 
fi = 20 Oe. 



Fig. 7 shows examples of vortex configurations ob- 
served for vorticities L = 2 to 13. The images shown 
in Fig. 7 were obtained in several different experiments 
and on samples of different sizes (see figure caption). We 
note that the same vorticity L could be obtained for dif- 
ferent combinations of the applied field and the size of the 
square, e.g., i^ = 6 was found for if = 60 Oe, a = 2 jam 
and if = 40 Oe, a = 2.5 jam - see images in Figs. 8(b) 
and 7(e), respectively. Sometimes two different vortic- 
ities were found in the same experiment for nominally 
identical squares, e.g., both i^ = 9 and i^ = 10 were 
found for if = 35 Oe and a = 3.5 /xm - see images in 
Figs. 7(g), (h), (i). The latter finding can be explained 



by slightly different shapes of individual squares or by 
an extra vortex captured during field cooling - see Ref.^« 
for a more detailed discussion, where the same effect was 
found for circular mesoscopic disks. Overall, the vorticity 
as a function of the applied field if showed the same be- 
haviour as that found earlier for circular disks^, i.e., the 
square dots showed strong diamagnetic response for small 
vorticities L < 10 (also observed earlier in disks with a 
strong disorder—) while for larger vorticities the extra 
demagnetisation per vortex saturated at d^/^ ~ 0.2, in 
excellent agreement with earlier numerical studies^^. 

Most of the vortex configurations shown in Fig. 7 rep- 
resent just one of several possible states for each vor- 
ticity (with the exception of images (h) and (i) which 
both correspond to L = 10). Indeed, for most vortici- 
ties we found more than one well-defined vortex config- 
uration and some of these were found with almost the 
same probability, indicating that, in agreement with the- 
ory described above, vortices in mesoscopic squares form 
not only the ground, but also metastable states, and the 
energies of the latter are often very close to the energy of 
the ground state. This conclusion follows from our statis- 
tical analysis of all observed vortex configurations which 
resulted in histograms such as those shown in Fig. 8 for 
L = 2, 4, 5, and 6. For L = 2 and 4, the most frequently 
observed states agree with the ground states found the- 
oretically (see Fig. 3(b), (d)) and the metastable states 
appear to have similar energies, as they are found with 
similar probabilities. As expected, both states for L = 2 
and two of the states for L = 4 have vortices sitting 
along the symmetry axes of the square, with the diag- 
onal axis being slightly preferable. The third state for 
L = 4 (on the right-hand side in Fig. 8(a)) is more un- 
usual in that the vortices are sitting in the apexes of 
a rhombus that is slightly rotated with respect to the 
diagonal of the square. Although this particular state 
did not come out in the numerical simulations^^, it was 
found with a high probability in experiment and, more- 
over, the rhombus-based vortex configurations were also 
found for larger vorticities both in experiment (see, e.g.. 
Fig. 7(1) for L = 13) and theory (see rhombic inner shells 
for i^ = 12 and 16 in Figs. 3(1) and 4(b), respectively). 

For L = 6, one of the two most frequently observed 
states (also shown in Fig. 7(e)) corresponds exactly to 
the ground state found numerically (Fig. 3(f)) but the 
state found in experiment with the highest probability 
is the more symmetric two-shell configuration with the 
outer shell having the same pentagon shape as that found 
for i^ = 5. This L = 6 state can be viewed as a direct 
precursor of the two-shell states for i^ = 7 and 9, which 
were found as ground states both in theory (Fig. 3(g),(i)) 
and experiment (Fig. 7(f), (g)). For L = two possi- 
ble states - a two-shell configuration with one vortex in 
the center (1,4) and four vortices in the corners and a 
pentagon-like configuration (5) - were found in experi- 
ment and in numerical simulations. However, numerical 
simulations found a slightly lower energy for the two-shell 
configuration (1,4) (see Fig. 5), while in experiment the 
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pentagon-shaped configuration was found to appear more 
frequently. Tiiis discrepancy is unlikely to be related to 
the non-ideal character of the experimental squares: As 
we show below, neither the roughness of the boundaries, 
nor the presence of some pinning in the experimental 
samples have any noticeable effect on the observed vor- 
tex configurations, due to strong confinement (see, e.g.. 
Fig. 2). It is possible that, due to the very small dif- 
ference in free energies between the two states (which 
becomes practically negligible for samples with a ^ A), 
the vortex configurations for L = 5 are particularly sen- 
sitive to the exact sample size (in experiment the squares 
are almost 10 times larger than in the analysis of Fig. 5). 
The sensitivity of vortex configurations to sample size 
was studied in detail for circular disks (see Ref. 9) and 
was indeed found to affect the stability of some (but not 
all) vortex states. For higher vorticilties, L = 7 to 13, 
we found well defined two-shell configurations most of 
which correspond to the stable configurations found nu- 
merically. The outer shell in these configurations was 
either square (see Figs. 7(g)- (k) for L = 9 to 12), circular 
(L = 7, Fig. 7(f)) or rhombic {L = 13, Fig. 7(1)) with 
vortices of the inner shell either sitting along one of the 
symmetry axes of the square, as for L = 2, or forming a 
triangle, as for L = 3. For certain matching vorticities 
(L = 9 and 12), the observed two-shell configurations 
correspond to a square vortex lattice. 

We note that the irregularities of the sample shape and 
uneven boundaries of some of our dots have, surprisingly, 
no discernable effect on the observed configurations of 
vortices (i.e., the vortices form regular, symmetric pat- 
terns). For example, the dots in Figs. 7(j),(k) have espe- 
cially rounded corners and very rough boundaries but the 
vortex configurations have square symmetries. Similarly, 
the same L = 6 state was found in dots with rounded 
corners, as in Fig. 7(e), and in almost perfect squares, as 
in the image shown in Fig. 8(b). Furthermore, we found 
that for a given value of L the observed configurations 
did not depend on the sample size or the applied field, 
at least within the studied field range see Fig. 9 for an 
example. 

Finally, we compared the experimentally observed po- 
sitions of vortices within the square dots with those 
found numerically and found an excellent agreement, as 
demonstrated by Fig. 9. Here we show a superposition 
of theoretical images from Fig. 3 and experimental im- 
ages for the same vortex configurations. Two of the im- 
ages (Figs. 9(d),(e)) compare the same theoretical con- 
figuration with experimental images obtained on dots of 
different sizes in different applied fields {H = 40 Oe, 
a = 2.5 /im and = 60 Oe and a = 2 /im, respectively) 
illustrating the point made above that the vortex config- 
urations do not depend on the sample size and/or applied 
field. 

Overall, despite the inevitable presence of some disor- 
der in our samples, which was not taken into account in 
the calculations, there is a very good agreement between 
the observed vortex configurations and the calculated 




L=2 L=4 




L=5 L=Q 



FIG. 8: Histograms of different vortex states observed for 
vortcities L = 2, 4 (for squares with a = 2/jm) (a) and 
L = 5 (for squares with a — 2/jm) and 6 (b) (a — 2/jm 
and a = 2.5/im). SEM images of the corresponding vortex 
configurations are shown as insets. 



vortex patterns. The main features of the vortex states 
revealed by experiment is formation of vortex shells with 
predominantly square symmetry for vorticities L > 7 and 
vortex patterns following the main symmetry axes of the 
square for small vorticities L < 4. The two intermedi- 
ate vorticities L = 5 and 6 appear to be a special case: 
Here the mismatch between the square shape of the dot 
and the natural symmetry of the vortex lattice is more 
difficult to accommodate and the preferred vortex con- 
figurations turned out to be the pentagon-shaped shell 
for L = 5 and three different patterns for L = 6, none of 
which has the four-fold symmetry of the square. 



V. CONCLUSIONS 

We performed a systematic study of vortex configu- 
rations in mesoscopic superconducting squares and com- 
pared the results with vortex patterns observed experi- 
mentally in /im-sized Nb squares using the Bitter deco- 
ration technique. 
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FIG. 9: (Color online) Comparison of the experimentally ob- 
served positions of vortices within the square dots with those 
found numerically. Superimposed on the experimental images 
are vortex configurations shown in Figs. 3(c),(d),(e),(f),(g). 
Two experimental images for L = 6 ((d) and (e)) are super- 
imposed on the same theoretical image (Fig. 3(f)), to demon- 
strate that the observed configurations did not depend on the 
sample size or the applied field (for image (d) H = 40 Oe, 
a ~ 2.5 /xm, for image (e) H = 60 Oe, a ~ 2 /xm). 



In the theoretical analysis we relied upon the analytical 
solution of the London equation in mesoscopic squares by 
using the Green's function method and the image tech- 
nique. The stable vortex configurations were calculated 
using the technique of molecular-dynamics simulations 
simulating the stimulated annealing process in experi- 



ments. 

We revealed the filling rules for squares with growing 
number of vortices L when gradually increasing the ap- 
plied magnetic field. In particular, we found that for 
small L vortices tend to form patterns that are commen- 
surate with the symmetry of the square boundaries of 
the sample. The filling of "shells" (similar to mesoscopic 
disks) occurs by periodic filling of the outermost and in- 
ternal shells. With increasing vorticity, the outermost 
shell is filled until it is complete (i.e., the number of vor- 
tices in it becomes 4A^, where N is an integer, i.e., com- 
mensurate with the square boundary). Then vortices fill 
internal shells untill the number of vortices becomes large 
enough to create the outermost shell with 4:{N + 1) vor- 
tices. Again, after that vortices fill internal shells. With 
increasing vorticity, the shell structure becomes less pro- 
nounced, and for large enough L the vortex patterns in 
squares becomes a traingular lattice distorted near the 
boundaries. 
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